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ABSTRACT 

Differential rotation in stars generates toroidal magnetic fields whenever an initial seed poloidal field is 
present. The resulting magnetic stresses, along with viscosity, drive the star toward uniform rotation. 
This magnetic braking has important dynamical consequences in many astrophysical contexts. For 
example, merging binary neutron stars can form "hypermassive" remnants supported against collapse 
by differential rotation. The removal of this support by magnetic braking induces radial fluid motion, 
which can lead to delayed collapse of the remnant to a black hole. We explore the effects of magnetic 
braking and viscosity on the structure of a differentially rotating, compressible star, generalizing our 
earlier calculations for incompressible configurations. The star is idealized as a differentially rotating, 
infinite cylinder supported initially by a polytropic equation of state. The gas is assumed to be 
infinitely conducting and our calculations are performed in Newtonian gravitation. Though highly 
idealized, our model allows for the incorporation of magnetic fields, viscosity, compressibility, and 
shocks with minimal computational resources in a 1+1 dimensional Lagrangian MHD code. Our 
evolution calculations show that magnetic braking can lead to significant structural changes in a star, 
including quasistatic contraction of the core and ejection of matter in the outermost regions to form 
a wind or an ambient disk. These calculations serve as a prelude and a guide to more realistic MHD 
simulations in full 3+1 general relativity. 

Subject headings: gravitational waves — MHD — stars: neutron — stars: rotation 



1. INTRODUCTION 

Magnetic fields play a crucial role in determining the 
evolution of many relativistic objects. Some of these 
systems are promising sources of gravitational radiation 
for detection by laser interferometers now under design 
and construction, such as LIGO, VIRGO, TAMA, GEO, 
and LISA. For example, merging neutron stars can form 
a differentially rotating "hypermassive" remnant (Rasio 
& Shapiro 1992, 1999; Baumgarte, Shapiro, & Shibata 
2000; Shibata & Uryu 2000). These configurations have 
rest masses that exceed both the maximum rest mass 
of nonrotating spherical stars (the TOV limit) and uni- 
formly rotating stars at the mass-shedding limit ( "supra- 
massive stars"), all with the same polytropic index. This 
is possible because differentially rotating neutron stars 
can support significantly more rest mass than their non- 
rotating or uniformly rotating counterparts. Baumgarte, 
Shapiro & Shibata (2000), have performed dynamical 
simulations in full general relativity to demonstrate that 
hypermassive stars can be constructed that are dynam- 
ically stable against radial collapse and nonradial bar 
formation. The dynamical stabilization of a hypermas- 
sive remnant by differential rotation may lead to delayed 
collapse and a delayed gravitational wave burst. The rea- 
son is that the stabilization due to differential rotation, 
although expected to last for many dynamical timescales 
(i.e. many milliseconds), will ultimately be destroyed 
by magnetic braking and/or viscosity. These processes 
drive the star to uniform rotation, which cannot support 
the excess mass, and will lead to catastrophic collapse, 
possibly accompanied by some mass loss. 
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Baumgarte & Shapiro (2003) discuss several other rel- 
ativistic astrophysical scenarios in which magnetic fields 
will be important. Neutron stars formed through core 
collapse supernovae are probably characterized by sig- 
nificant differential rotation (see, e.g., Zwerger & Miiller 
1997; Rampp, Miiller, & Ruffert 1998; Liu & Lindblom 
2001; Liu 2002, and references therein). Conservation 
of angular momentum during the collapse is expected 
to result in a large value of (3 = T/\W\, where T is 
the rotational kinetic energy and W is the gravitational 
potential energy. However, uniform rotation can only 
support small values of (3 without leading to mass shed- 
ding (Shapiro & Teukolsky, 1983). Thus, nascent neu- 
tron stars from supernovae probably rotate differentially, 
giving magnetic braking and viscosity an important dy- 
namical role. Short period gamma-ray bursts (GRB's) 
are thought to result from binary neutron star merg- 
ers (Narayan, Paczynski, & Piran 1992) or tidal disrup- 
tions of neutron stars by black holes (Ruffert & Janka 
1999). Meanwhile, long period GRB's are believed to 
result from collapses of rotating massive stars to form 
black holes (MacFadyen & Woosley 1999). In most cur- 
rent models of GRB's, the burst is powered by rotational 
energy extracted from the neutron star, black hole, or 
surrounding disk (Vlahakis & Konigl 2001). This en- 
ergy extraction can be accomplished by strong magnetic 
fields, which also provide an explanation for the collima- 
tion of GRB outflows into jets (Meszaros & Rees 1997; 
Sari, Piran, & Halpern 1999; Piran 2002) and for the ob- 
served gamma-ray polarization (Coburn & Boggs 2003). 
The dynamics of supermassive stars (SMS's), which may 
have been present in the early universe, will also depend 
on magnetic fields if the SMS's rotate differentially. Loss 
of differential rotation support will affect the collapse of 
SMS's, which has been proposed as a formation mech- 
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anism for supermassive black holes observed in galaxies 
and quasars, or their seeds (see Rees 1984, Baumgarte 
& Shapiro 1999, Bromm & Loeb 2003, and Shapiro 2003 
for discussion and references). Finally, the effectiveness 
of the r-mode instability in rotating neutron stars may 
depend on magnetic fields. This instability has been pro- 
posed as a mechanism for limiting the angular velocities 
of neutron stars and for producing quasi-periodic grav- 
itational waves (Andersson 1998; Friedman & Morsink 
1998; Andersson, Kokkotas, & Stergioulas 1999; Lind- 
blom, Owen, & Morsink 1998). Due to flux- freezing, 
however, the magnetic fields may distort or suppress the 
r-modes (see Rezzolla et al. 2000, 2001a, b and refer- 
ences therein). 

Motivated in part by the growing list of important, un- 
solved problems which involve hydromagnetic effects in 
strong-field dynamical spacetimes, Shapiro (2000, here- 
after Paper I) performed a simple, Newtonian, MHD cal- 
culation to highlight the importance of magnetic fields 
in differentially rotating stars. In this model, the star is 
idealized as a differentially rotating infinite cylinder con- 
sisting of a homogeneous, incompressible conducting gas 
in hydrostatic equilibrium. The magnetic field is taken 
to be purely radial initially and is allowed to evolve ac- 
cording to the ideal MHD equations. A constant shear 
viscosity is also allowed. Note that, though highly ide- 
alized, this model allows several important physical ef- 
fects to be included in a 1+1 evolution code, requiring 
much less computational effort than more realistic 3+1 
simulations. 

Paper I presents results for the behavior of differen- 
tially rotating cylinders in three physical situations: zero 
viscosity with a vacuum exterior, nonzero viscosity with 
a vacuum exterior, and zero viscosity with a tenuous am- 
bient plasma atmosphere. In the first case, for which the 
solution is analytic, the magnetic field gains a toroidal 
component which oscillates back and forth in a standing 
Alfven wave pattern. The angular velocity profile os- 
cillates around a state of uniform rotation, with uniform 
rotation occurring at times when the azimuthal magnetic 
field is at its maximum magnitude. At these times, a sig- 
nificant amount of the rotational energy has been con- 
verted to magnetic energy. In the absence of a dissipa- 
tion mechanism, these oscillations continue indefinitely. 
In the second case (vacuum exterior with nonzero vis- 
cosity), the oscillations are damped and the cylinder is 
driven to a permanent state of uniform rotation with a 
significant fraction of its initial rotational energy hav- 
ing been converted into magnetic energy and finally into 
heat. In the final case (zero viscosity with a plasma at- 
mosphere), Alfven waves are partially transmitted at the 
surface and carry away significant fractions of the an- 
gular momentum and energy of the cylinder. In all of 
these situations, only the timescale, and not the quali- 
tative outcome, of the dynamical behavior depends on 
the strength of the initial magnetic field. Thus, even for 
a small initial seed field, magnetic braking and viscosity 
eventually drive the cylinders toward uniform rotation. 

The braking of differential rotation by magnetic fields 
is also currently being studied in a spherical, incompress- 
ible neutron star model (Liu & Shapiro, in preparation). 
This calculation treats the slow-rotation, weak-magnetic 
field case in which T <C i? mag <C |W|, where E mag is the 
total magnetic energy. Consequently, the star is spheri- 



cal to a good approximation and poloidal velocities in the 
9~ and r— directions can be neglected on the timescales 
of interest. Liu & Shapiro solve the MHD equations for 
the coupled evolution of the angular velocity fl and the 
azimuthal magnetic field in both Newtonian gravity 
and general relativity. The resulting angular velocity and 
poloidal field profiles often develop rich small-scale struc- 
ture due to the fact that the eigenfrequencies of B^ vary 
with location inside the star. 

In many cases, the loss of differential rotation support 
due to magnetic braking is expected to lead to interest- 
ing dynamical behavior. In hypermassive stars, magnetic 
braking may lead to catastrophic gravitational collapse. 
In general, however, one expects a range of radial mo- 
tions, including oscillations, quasistatic contraction, and 
ejection of winds. These behaviors are not present in 
the results of Paper I because of the assumptions of ho- 
mogeneity and incompressibility. In that case, no radial 
velocities occur as the system evolves away from equi- 
librium. In the present paper, we explore the effects of 
compressibility by generalizing the calculations of Paper 
I to Newtonian cylindrical polytropes. Though other as- 
pects of the present model are again highly idealized, 
allowing for compressible matter results in a much more 
realistic spectrum of dynamical behavior, including ra- 
dial contraction, mass ejection, and shocks. By formu- 
lating the problem in 1+1 dimensions and working in 
Lagrangian coordinates, we are able to solve the MHD 
equations essentially to arbitrary precision. Our simula- 
tions again serve to identify most of the key physical and 
numerical parameters, scaling behavior, and competing 
timescales associated with magnetic braking and differ- 
ential rotation. The structure of this paper is as follows: 
SectionEldescribes our basic model. In Section|3| we dis- 
cuss the structure of equilibrium cylindrical polytropes 
which serve as initial data for our evolution calculations. 
In Section^] we set out the fundamental MHD evolution 
equations and put them in a convenient nondimensional 
form. Section [5] describes the results of MHD evolution 
calculations for several choices of the parameters which 
serve to illustrate the effects of compressibility and dif- 
ferential rotation. We discuss the significance of these 
calculations in Sectional 

2. BASIC MODEL 

As discussed in Paper I, differential rotation in a 
spinning star twists up a frozen-in, poloidal magnetic 
field, creating a very strong toroidal field. This pro- 
cess will generate Alfven waves, which can redistribute 
the angular momentum of a star on timescales less than 
~ 100 s (10 12 G/i?o), where Bo is the characteristic initial 
poloidal field. Shear viscosity also redistributes angu- 
lar momentum. However, molecular viscosity in neutron 
star matter operates on a typical timescale of ~ 10 9 s in a 
star of ~ 10 9 K, so it alone is much less effective in bring- 
ing the star into uniform rotation than magnetic brak- 
ing, unless the initial magnetic field is particularly weak. 
Turbulent viscosity, if it arises, can act more quickly. 

We wish to track the evolution of differential rotation 
and follow the competition between magnetic braking 
and viscous damping. Following Paper I, we model a 
spinning star as an infinite, axisymmetric cylinder with 
a vacuum exterior. We adopt cylindrical coordinates 
(r, (/>, z), with the z— axis aligned with the rotation axis 
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of the star, and assume translational symmetry in the z— 
direction. We take the magnetic field to have components 
only in the r— and directions. The fluid is assumed 
to be perfectly conducting everywhere (the ideal MHD 
regime). In some of our models, we allow for the pres- 
ence of a constant shear viscosity. For protoneutron stars 
with T > lO 9 !^ , the coefficient of bulk viscosity, due to 
the time lag for reestablishing beta equilibrium following 
a change in density, becomes comparable to the shear 
viscosity coefficient (Sawyer 1989). We do not account 
for bulk viscosity, however, since our models are strongly 
differentially rotating (and hence have large amounts of 
shear), while large scale density changes in our evolution 
calculations occur quasistatically (see Section [SJ. 

To construct initial data for our dynamical simulations, 
we assume a polytropic equation of state (EOS): P = 
Kp T . The equations of hydrostatic equilibrium are then 
solved with an assumed form for the initial differential 
rotation law (discussed in Section [5}. We take the <f>— 
component of the magnetic field to be zero initially. We 
then evolve the system by numerically integrating the full 
set of coupled MHD equations. We allow for heating due 
to shocks and viscous dissipation during the evolution, 
and assume a T— law EOS, P = {T — l)ps, where s is the 
specific internal energy. 

3. EQUILIBRIUM CYLINDRICAL POLYTROPES 

Here we discuss the initial data for our simulations. We 
derive the equations for the static case, and then general- 
ize to include rotation and radial magnetic fields. Finally, 
we derive a virial relation for equilibrium infinite cylin- 
drical stars, discuss the stability of these models to radial 
perturbations, and construct some numerical models. 

3.1. Static Polytropes 

The fundamental equations for a static polytrope are 
the polytropic EOS 

P = Kp r = Kp 1+ i, (1) 

where p is the density, P is the pressure, K and T are 
constants, and n is the polytropic index, and the equation 
of hydrostatic equilibrium 

-VP = -V$, (2) 
P 

where $ is the Newtonian gravitational potential, which 
satisfies Poisson's equation, 

V 2 $ = Airp. (3) 

Note that we will work in units such that G = 1 . To find 
V$, we first define p(r) to be the mass per unit length 
interior to radius r. This gives 



(4) 



p,(r) = 2tt / p{r )r dr 



Integrating equation @ and imposing regularity at the 
origin yields 

<*g _ w (5) 

dr r 

Differentiating equation J5J in cylindrical radius r and 
combining with equation @ yields 



d (rdP\ 

Tr lp* ' = ~ Anp{r)r - 



(0) 



This equation can be generalized for infinite planar or 
spherical geometries through similar arguments. The re- 
sult is 

„a-l jp^ 

' = -4np(r)r a - L , (7) 



d_ 
dr 



P 



dr 



where a = 1,2,3 corresponds to planar, cylindrical, and 
spherical geometry respectively. 

We will make use of a common nondimensionalization 
of equation Q) to facilitate numerical solutions. We de- 
fine nondimensional quantities according to 



P = p c o 
r = a£ 



(8) 



(n + l)Kp£ /r 
An 



1/2 



where p c is the central density. Equation (JJJl then be- 
comes a generalized Lane-Emden equation, 



1 



e*- 1 dn 



(9) 



Hereafter we will treat a — 2 (cylindrical geometry). 
Spherical polytropes were described extensively by Chan- 
drasekhar (1939); planar polytropes were constructed by 
Shapiro (1980). The boundary conditions on equation 
Q are easily obtained from physical arguments. From 
equation (JSJ, it is evident that 9(0) = 1. The condition 
on 9' is obtained by using the fact that p — * p c so that 



equation Q implies p(r) cx 



0, hence equations 



d), ©, and © yield dp/dr = 0, and hence 0'(O) = 0. 
The equations should be integrated from the origin out 
until 9 = 0, which implies P = 0. This point is the sur- 
face of the star and will be denoted by £i. This gives 
the pressure and density profiles for equilibrium configu- 
rations. Some analytic solutions to equation @ are, for 



0. 



1 - £ 2 /4 and for n = 1, = Jo(£) ; where Jo 



is the Bessel function of order zero. 

3.2. Rotating Cylindrical Polytropes 

We will now extend the cylindrical polytrope model 
to include the effects of rotation. The hydrostatic equi- 
librium equation is modified by a term accounting for 
centrifugal acceleration. The result is 

1 dP d$ 2 

+ n 2 (r)r, (10) 



p dr 



dr 



where £l(r) is the angular velocity at radius r (the in- 
dependence allows for differential rotation). Performing 
the same manipulations used to get equation © gives 



d frdP\ 



dr \p dr J 



= 2VL z r + r z - 



,dn 2 

dr 



Anpr. (11) 



Using the same nondimensionalization as equation |JS} 
and making the definition 



n 2 

KPc 



(12) 



gives the cylindrical Lane-Emden equation with rotation, 



4d£ 



0. 



(13) 
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By identical reasoning, this equation has the same 
boundary conditions as equation 0. 

There are two analytic solutions for uniform rotation 
(( = constant). The analytic solutions are for n = 0, 
0(f) = 1 - (1 - C/2)£ 2 /4 and for n = 1, 0(f) = (1 - 
C/2) Jo(0 + C/2- In the rest of the paper, we will adopt 
the following profile 



n(r) = 



fin 



1 



cos 



Tvr 

"R2 



2 \ 1 



(14) 



where R is the radius of the star. Letting ( = Vt^jnp,. 
gives the nondimcnsional form 



C(f) 



4 



1 



(15) 



Though this rotation law was chosen arbitrarily, it sat- 
isfies the boundary condition given below (eq. |45j ) 
and has the expected property that f2 decreases mono- 
tonically with r. Shibata and Uryu (2000) found that 
hypermassive binary merger remnants typically have 
Q(0)/Q,(R) 3. Slightly modifying our rotation law to 
satisfy this condition resulted in a set of models qual- 
itatively similar to those discussed below, usually with 
slightly larger allowed values of j3 = T/\W\. To solve 
equation l|13|l for rotation profile of equation l)15[l. we 
must iterate our solutions to equation (fH5)) to identify 
fi, the surface radius of our star. As will be shown later, 
these equilibrium solutions are also valid for cylinders 
with radial magnetic fields, so we will use these solu- 
tions as initial data for our dynamical evolutions. We 
remark how simple it is to construct differentially rotat- 
ing, infinite cylindrical polytropes, described by ordinary 
differential equations in one dimension, in comparison to 
finite rotating stars, which are defined by partial differen- 
tial equations in two dimensions (see e.g. Tassoul 1978). 

3.3. Virial Theorem and Stability 

We will use the virial theorem to check the results 
of our equilibrium calculations and, later on, to test 
whether a configuration has reached a new equilibrium 
state following dynamical evolution. With the defini- 
tions 



W = - px- V$d 3 x, 



T=~ I pv'd'x, 

n= / Pd 3 x , 



(16) 



the virial theorem for an infinite cylindrical equilibrium 
star is 

2T + VK + 2II = 0. (17) 

Notice that the numerical coefficients in this equation dif- 
fer from the analogue for bounded configurations. (For 
a derivation, see Appendix El) Also note that as shown 
in Appendix even if magnetic fields are present, this 
equation still holds unmodified for equilibrium stars. In 
our case it is convenient to write d 3 x = LdA — 2irLrdr, 
where L is an arbitrary length in the z-direction. We 
divide each quantity in equation l|l(jf> by L to obtain 
quantities per unit length (we will continue to use the 



symbols W, T, and II). Hereafter, we will refer to all 
extensive quantities as quantities per unit length in the 
z— direction. These definitions give an identical virial 
theorem to equation i|17|) . Using equations J2J| and @, 
the definition of W can be integrated analytically to give 



W 



(18) 



where p, t — fJ-(R) is the total mass per unit length of 
the star. 

The virial theorem can be used to find an upper limit 
on P = T/\W\ for any equilibrium cylindrical star. Since 
the pressure is always positive, II is positive definite. 
Similarly, p? t is positive definite, so W is negative or 
W = — \W\. Using these facts in equation (|T7jl implies 
that for any equilibrium star 



1 



(19) 



We require our equilibrium solutions to be initially sta- 
ble to radial perturbations to ensure that we are see- 
ing the secular effects of magnetic fields in the evolution 
and not a hydrodynamical radial instability. For cylin- 
ders, the critical Y for stability is Y > 1, permitting any 
n > (sec Appcndix[EJ). For spheres, however, one needs 
r > 4/3, which in turn requires n < 3. Note that rota- 
tion lowers the critical Y, which does not change the limit 
on n, since any realistic EOS allows us to restrict our- 
selves to n > 0. Also note that the presence of a radial 
magnetic field does not affect this stability analysis, as 
shown in Appendix 151 

3.4. Equilibrium Models 

Here we present the results of our equilibrium model 
calculations. We write equation ljl3[) as a system of two 
coupled first order ODE's and integrate using a fourth 
order Runge-Kutta method. Below we discuss the com- 
parison with an analytic Roche model for centrally con- 
densed rotating cylinders and give a table of data for 
cylindrical polytropes. 

3.4.1. Roche Model 

The Roche model holds for uniformly rotating, cen- 
trally condensed, polytropic cylinders. Our discussion 
is patterned after Shapiro & Teukolsky (1983). We will 
show below that this approximate model improves with 
increasing n, as expected. Equation (|10|l may be writ- 
ten as 

-WP + V($ + $ c ) = 0, (20) 
P 

where $ c = — £l 2 r 2 /2. Because of the central condensa- 
tion, near the surface fj, — ► /i t and 



9$ 2u t ^ „ , 
or r 



The specific enthalpy of a polytrope is given by 
, . 1 dP Y P 

ir ~in = 

p or 



h(r) 



(21) 



(22) 



Y-lp 

Equation 120fl can now be integrated to obtain 

h(r) + $(r) + $ c (r) - K, (23) 

where K is an integration constant. Since the star is very 
centrally condensed and equation ll'.'il) holds throughout 
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Table 1 

Static Cylindrical Polytropes 



n 


Ci 


(f k 


pc/p 


0.0 


2.0 


2.0 


1.0 


1.0 


2.405 


1.248 


2.316 


2.0 


2.921 


0.925 


4.611 


3.0 


3.574 


0.740 


8.629 


5.0 


5.428 


0.532 


27.67 



the envelope, K is assumed to be the same as in the 
nonrotating case. Let Ro be the radius of the nonrotating 
star. Then, since the enthalpy vanishes at the surface, 



$(JZo) = K= 2 A i t lni? . 



(24) 



A stable rotating star (i.e. not shedding mass) must 
have a radius R\ such that h(R%) = 0. Let the effec- 
tive potential be defined as <E> e ff = $ + and let the 
maximum of the effective potential be $ mal . Then us- 
ing equation i(23() and the fact that h(r) has a zero, this 
stability condition becomes 

$max > K, (25) 

where <l> max is determined by setting d$ c s/dr = 0, giving 



Inserting this value into $ n 
ity condition 

- Ht + 2^ t 



(26) 

gives the rotational stabil- 



^ 2/it In Rq . 



(27) 



This implies that to avoid mass shedding, the star 
must satisfy 

(28) 



^ eR? 



Then in terms of nondimensional variables of the nonro- 
tating star, the stability criterion in the Roche approxi- 
mation becomes 



1 4 89 

C s£ 



(29) 



£=6 



employing equations <|I 2(1 . and l)28[l. Now, we use 
our code to find the maximum £ for solutions to equa- 
tion 113(1 and compare to the Roche model results. We 
denote the maximum £ from our code as Q m and the 
maximum £ from the Roche model, equation 1(29(1 . as 
Ci?. We find the following sets (n, C m , ( R ): (1.0, 0.57425, 
0.31767), (3.0, 0.10948, 0.08527), (5.0, 0.02977, 0.02659), 
and (10.0, 0.001795, 0.001767). These results confirm 
that the Roche model becomes increasingly more accu- 
rate as n increases and the degree of central concentra- 
tion increases. 

3.4.2. Numerical Results 

Table 1 summarizes numerical results for static cylin- 
drical polytropes. As n increases, the star becomes more 
centrally condensed. However, this is much less pro- 
nounced than in the spherical case, where n = 3 gives 
Pc/p — 54.18 compared to p c /p — 8.629 for cylinders. 



0.5 



0.4 



0.3 







0.2 



0.1 



IV 



n 

Fig. 1. — Maximum (3 = T/\W\ attainable for uniformly rotat- 
ing stars (mass shedding limit; the solid line) and the values of f) 
which we will evolve numerically (marked with crosses and labeled 
I- VI). Note that the virial theorem sets the limit [3 < 0.5 for any 
equilibrium configuration. 



Table 2 gives the key results for uniformly rotating cylin- 
drical polytropes spinning at the mass-shedding limit. 
Note that we have not found a maximum n for solutions 
with finite p t in the case of cylindrical polytropes. In 
contrast, polytropic spheres with n > 5 extend to infi- 
nite radius with infinite mass (Chandrasekhar 1939). An 
important consequence of our differential rotation law is 
that for n > 1.2, we can get a larger value of (3 than 
the maximum [3 for uniform rotation, /3 max . Since we ex- 
pect viscosity and magnetic braking to drive differential 
rotation toward uniform rotation, we expect the most 
interesting cases to be those where the star is forbidden 
to simply relax to uniform rotation. Figure ^ shows a 
plot /3 max versus n, where we have marked the cases we 
dynamically evolve below after incorporating magnetic 
fields and viscosity. We expect the cases falling below 
the curve to relax to a uniformly rotating equilibrium 
polytrope and the cases above the curve to undergo ma- 
jor changes away from polytropic behavior. 

4. MHD EQUATIONS 

We now assemble the fundamental equations for our 
MHD simulations. First, the fluid obeys the continu- 
ity equation, 

^+V-(pv) = 0. (30) 

The fluid motion is governed by the magnetic Navier- 
Stokes equation: 



<9v 

w + (v.V)v : 



1 

-VP 

P 
V 



(V x B) x B 



4-7T/9 



(V 2 v+ IV(V-v)), 



(31) 



where B is the magnetic field and r\ is the coefficient of 
dynamic viscosity, which is related to the familiar kine- 
matic viscosity according to v = r//p. We note that if 
B = — rj, and v = £1 x r, fi = ilz, equation 1(31(1 



G 



Table 2 

Maximally Rotating Cylindrical Polytropes 
Uniform Rotation 3 Differential Rotation 13 



n 


0max 


^rnax / yj ^ Pc 


0max 


^max j Pc 


0.0 


0.500 


1.414 


0.086 


1.414 


0.01 


0.496 


1.402 


0.087 


1.414 


1.0 


0.228 


0.758 


0.179 


1.414 


3.0 


0.062 


0.331 


0.299 


1.414 


5.0 


0.015 


0.172 


0.368 


0.382 



a The maximum corresponds to the mass shedding limit. 

b Assumes the particular rotation law Q = ^p- ^1 + cos ^-^-^1 and a mono- 
tonically decreasing density profile. 



reduces to the hydrostatic equilibrium equation for ro- 
tating cylinders (eq. ^H])- The magnetic field B satisfies 
Maxwell's constraint equation 

VB = 0, (32) 
as well as the flux-freezing equation 
r)B 

— = V x (v x B). (33) 

Solving equation for the magnetic field together 
with the flux-freezing condition l|33[) requires that the 
radial component of the magnetic field be independent 
of time and given by 

B r = t > 0, (34) 

r 

where Bo is the value of the field at the stellar surface 
at t = 0. Although the magnetic field given by equation 
l|34(l exhibits a static line singularity along the axis at 
r = 0, it does not drive singular behavior in the fluid 
velocity or nonradial magnetic field. As these quantities, 
which are the main focus here, remain finite and evolve in 
a physically reasonable fashion, the line singularity poses 
no difficulty and requires no special treatment. Also, 
since the magnetic term in equation (eq. |31j ) is propor- 
tional to the curl of B, the radial magnetic field does not 
enter the hydrostatic equilibrium equation. Thus, the 
presence of an initially radial magnetic field does not af- 
fect the initial data equilibrium models described in Sec- 
tion |3J The virial theorem derived by taking a moment 
of equation (|31() is identical to equation (|17|l . which ap- 
plies in the B = case. (The magnetic field terms cancel 
as described in Appendix lA"l) Note that the virial theo- 
rem only holds for equilibrium configurations and is not 
expected to be satisfied during a dynamical evolution. 

We next consider the properties of the fluid itself. The 
equation of state, P = (T — l)pe, allows the entropy of a 
fluid element to increase if there is heating due to shocks 
or shear viscosity. The evolution of the specific internal 
energy is described by the First Law of Thermodynamics, 



-P- 



T- 



ds 



(35) 



ds 

dt dt \p ) dt 

where the Lagrangian derivative is defined as d/dt = 
d/dt + v • V and s is the specific entropy. The rate of 
energy generation due to viscosity is given by 2 
,e?s 1 , dvi 



T- 



dt 



p"* 3 Ox, 



(36) 



2 Equations 1361 and 1371 are expressed in Cartesian coordinates 
for simplicity. 



where <r^ is the viscous stress tensor (Landau 
1998, eq. 15.3) 



<1 



dvj_ 

dxi 



3lJ dxJ 



Lifshitz 



(37) 



and where summation over repeated indices is implied. 

We will now write our system of evolution equations in 
component form in Lagrangian coordinates (we work in 
an orthonormal basis). First, we define j as the specific 
angular momentum: j = rv^, where is the <p— compo- 
nent of the velocity. We also explicitly use the gradient 
of $ given in equation (JSJ. Then equations H3Uf) . (|31fl . 
O, (33, and become 



dp 
~dt 
dv r 
~~dt 



di 
dt 



(rv r ) 
2/y 



p d 
r dr 
1 dP 
p dr 

4?7 d 

3p dr 

4^™ 



3 

h 

r r 
1 d(rv 



d 



8npr 2 dr 



dBj, 
dt 


or 


de 




dp 


~di~ 


P 2 ( 






9p 


[(" 







rB, 



dr 
r\r d 
p dr 

dr 



ldj_ 

r dr 



(38) 

(39) 
(40) 
(41) 



rjr 
P 



, dv r 

dr 

1 ( d(rv r ) 
dr 



d_ 

dr 

2 



d_ 

dr 



(42) 



The terms on the r. h. s. in equation (|42|l involving deriva- 
tives of v r and j represent heating due to the presence 
of shear viscosity. We note that the heating terms which 
depend on v r would not be present for an incompressible 
fluid. To handle shocks, we supplement the pressure by 
an artificial viscosity term: P — > P + q, where the recipe 
for q is described in Appendix |0 

4.1. Initial Conditions and Boundary Values 

We assume that no azimuthal component of the mag- 
netic field is present initially, 

^(0,r)=0, r>0, (43) 
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recalling that B r is always given by equation (J33J. We 
are thus interested in the situation where the azimuthal 
field is created entirely by the differential rotation of the 
fluid, which bends the frozen-in, initially radial field lines 
in the azimuthal direction. Since the region outside of the 
cylinder is a vacuum, no azimuthal magnetic field can be 
carried into this region. Thus, must vanish at the 
surface: 

B (f> (t,R)=0, t>0. (44) 
Together with equation l|41|l . this implies 



0_ 

Or 



j(t,R) 



0, t>0. 



(45) 



The angular velocity — j/r 2 must be finite at the 
origin. This implies that j(t, 0) = for all time. Also, 
since the cylinder is axisymmetric, the radial velocity at 
the origin must vanish for all time: v r (t,0) = 0. Since 
we begin our evolutions with an equilibrium cylindrical 
polytrope model, we take the radial velocity to be ini- 
tially zero everywhere, that is, v r (0, r) — 0. 

We now consider the boundary conditions on the fluid 
variables at the free surface of the cylinder, r = R. Pres- 
sure and density must both vanish at this surface: 



p(t,R) = 
P(t,R) = Q, 



(46) 



t > 0. 



When a viscosity is present, the free surface imposes an 
additional condition. Let n be the unit normal vector to 
the surface. Then, at the free surface, balance of forces 
requires that 

- Pn t + a' i:j nj = 0, (47) 

(Landau & Lifshitz 1998, eq. 15.16). Since n is radial 
and P = on the surface of the cylinder, this condition 
becomes a' rr — a' r ^ = 0. Written in component form, 
a', — becomes equation (|4*5)l and a' rr = becomes 



d 1 

2—v r (t,R)--v r (t,R) 
or r 



0, t > 0. 



(48) 



4.2. Conserved Energy and Angular Momentum 

The magnetic Navier-Stokes equation (eq. j3J) admits 
two nontrivial integrals of the motion, one expressing 
conservation of energy and the other conservation of an- 
gular momentum of the star. Energy conservation may 
be expressed as 

E(t) = E kin (t) + E mag {t) + E lnt (t) + E gIav (t) 

= E(0), (49) 

where E k i n (t) is the kinetic energy of the matter, £ , mag (i) 
is the azimuthal magnetic energy, Ei nt (t) is the internal 
(thermal) energy, and E gTav (t) is the gravitational poten- 
tial energy. These contributions are as follows: 



£kin(*) = 5 J dAp{t,r) (v 2 r {t,r) 
E int (t)= I dA P (t,r)e(t,r), 



■^grav(^) ~~ 2 



dAp(t,r)$(t,r), 



(50) 



where $(t,r) is the gravitational potential, dA = 2irrdr, 
and all energies are taken per unit length as in Sec- 
tion Equation l|49|l is derived by rewriting the Eu- 
lerian time derivative of the total energy density using 
equations (|30|l . (|31|l . and lj^5(l and applying the diver- 
gence theorem. There is no explicit contribution from 
viscous heating in equation l|49|l because the internal en- 
ergy, £i n t, accounts for this. We note that E gTliv (t) is not 
the same as the gravitational term W appearing in the 
virial theorem for equilibrium cylinders, which is defined 
according to equation ifTBl) . By contrast, E gT!!lv (t) = W 
for spherical stars (Shapiro & Teukolsky 1983). Due to 
the cylindrical geometry of the present problem, W can- 
not be identified as the gravitational potential energy. 
However, the ratio T/|I-F|, involving quantities from the 
virial theorem, can still be taken as a measure of the 
degree of rotation for a given mass per unit length. 

Conservation of angular momentum per unit length is 
expressed as 



At) 



dAp(t,r)j(t,r) = J(0) 



(51) 



Equation H51J) is obtained by rewriting the Eulerian time 
derivative of the angular momentum density using equa- 
tions Q3U[) and l|31|l and applying the divergence theorem. 
We note that, consistent with the nonrelativistic MHD 
approximation, the electric field energy E 2 /8tt is not in- 
cluded in equation (|49|l and the angular momentum of 
the electromagnetic field /c 2 , where S is the Poynting 
vector, is not included in equation H51(l (Landau, Lifshitz 
& Pitaevskii 2000). 

The motivation for monitoring the conservation equa- 
tions during the evolution is twofold: physically, evalu- 
ating the individual terms enables us to track how the 
initial rotational energy and angular momentum in the 
fluid are transformed; computationally, monitoring how 
well the conservation equations are satisfied provides a 
check on the numerical integration scheme. 

4.3. Nondimensional Formulation 

We introduce nondimensional variables to facilitate 
comparison of the results of the various dynamical simu- 
lations we present in Sectional First we define an Alfven 
velocity using the scale of the radial magnetic field from 
equation (|34|l : 

- s 7§? <52) 

where, here and throughout, p c is the initial central den- 
sity. We define nondimensional quantities according to 

r = r/R, t = 2t/(R/v A ), f, = ir)/(p c v A R) (53) 
j = j/n Q R 2 , B^ = B+l ((4 m ) 1 / 2 (fi i?)) 

p = p/p c , v r = v r /2v A , i = e/{Anv 2 A ) 
A= p/{2np c R 2 ), P = P/(8irp c v 2 A ). 

Unless stated otherwise, we will work with nondimen- 
sional variables in all subsequent equations, but we omit 
the carets (") on the variables for simplicity. With the re- 
placement V = l/p and the relation {pr)^ 1 d / dr — d/dp, 
the system of evolution equations (eqs. EHJ-E3) ^> e ~ 
comes 

-ir = dp:^ ) (54) 



8 



dv r 
~~dt 



dt 

dt 

ds 
dt 



-2nr— - x— - — — (r 2 B 2 ) 
dn r 4r dfi * 



CoX 
2 



id_ 

2djl 



= -B 



j 2 r/r d ( 1 d(rv r ) 
r 3 6 dp, \r dr 
r/r 2 d ( dj 
4 dp \dr 2 
dv r 1 d ( j 



' dr 

P dp 
~i?dt 4 

1 2/7 

1 f d(rv r ) 
dr 



2 dr \ r 
CoX r/r 2 ( d 



(55) 
(56) 
(57) 



167T p 

dv r 
dr 



dr 



V r 

r 



8_ /Vr\ 

dr V r 2 / 



(58) 



Here, Co = C( r = 0) is as defined in equation ljT2T) . x = 
Slgi? 2 /2v A , and artificial viscosity is incorporated in P 
as described in Appendix IO 

5. DYNAMICAL EVOLUTION 

The MHD evolution cases which we will discuss in Sec- 
tions 1^21^31 ar e specified by four parameters: the poly- 
tropic index n, the rotation parameter Co = Q^/irp C) 
the Alfven speed va = B /^4:Trp c , and the coefficient 
of dynamic viscosity r/. We will give results for n — 
0.001, 1, 3, and 5 in order to probe the effects of the de- 
gree of central condensation. For n = 0.001 and n = 1, 
all differentially rotating cylinders have (3 smaller than 
the maximum value allowed for uniform rotation, /3 max . 
For these cases, Co is chosen so that a moderate value 
of P = T/\W\ is obtained. For differentially rotating 
cylinders with n = 3 and n = 5 however, it is possi- 
ble to choose values of Co f° r which (3 > /? max . Thus, 
for both of these cases, we give results for two values of 
/3, one with (3 > /3 max and one with (3 < /3 max . Real- 
istic cold neutron star EOS's are expected to be fairly 
stiff, corresponding to polytropes with 0.5 < n < 1.0. 
For polytropic spheres, this range of n corresponds to 

1.84 < p c /p < 3.29. To obtain the same range of p c /p 
for cylinders requires 0.70 < n < 1.49. Our n = 1 model 
lies in this range, but we will also study the more com- 
pressible models since these allow /3 > /3 max with our 
adopted angular velocity profile and since we wish to un- 
derstand the effects of compressibility. We also note that 
our results may be relevant for magnetic braking in radi- 
ation pressure dominated SMS's, which have n « 3 and 
Pc/p » 54. 

For each case below, va is chosen so that A4a/T = 

7.85 x 10~ 3 , where A4a gives the relative scale of the 
energy per unit length associated with the initial radial 
magnetic field: 



B 2 

M A = nR 2 ^ = \-Kp c R 2 v\. 



(59) 



Because of the static line singularity at r = 0, the ac- 
tual energy associated with the radial field B r is not well 
defined. For this reason, the energy scale A^a is used 
to characterize the radial field strength. We choose the 
ratio Ma/T to be the same for all of our runs because 



this ratio roughly determines the extent to which the 
radial field lines are wound up. Thus, keeping this ra- 
tio constant facilitates comparison between the different 
evolution cases. Choosing Ma *C 1 shows the conse- 
quences of an initially weak magnetic field in influencing 
the dynamical behavior of an equilibrium star. This is 
the situation of greatest astrophysical interest for core 
collapse, hypermassive neutron star evolution, and other 
relativistic MHD scenarios. We define an Alfven wave 
crossing time as £a = R/va, which in nondimensional 
units is just tA = 2. Thus, the basic time unit for the 
evolution code is the Alfven time. 

Lastly, for cases with shear viscosity, the coefficient of 
dynamical viscosity is chosen as r\ = 0.2 in nondimen- 
sional units. This intermediate value is chosen in order 
that the viscosity in the star be sufficiently large for dissi- 
pation to become evident in a few Alfven timcscalcs, but 
sufficiently small so that viscous damping of differential 
rotation does not completely suppress the growth of the 
azimuthal magnetic field. The characteristic viscous dis- 
sipation timescale can be taken as t v = 1 /rj in our nondi- 
mensional units, which is equivalent to t n = R 2 p c /8r) in 
physical units (see eq. [SI)- The choice 77 = 0.2 corre- 
sponds to the ratio of the viscous to Alfven timescales 
tn/tA = 2.5. 

To illustrate the hierarchy of timescales governing the 
evolution of relativistic rotating stars, we will evaluate 
the relevant timescales for neutron star parameters con- 
sistent with the stable hypermassive remnant found in 
the simulations of Shibata & Uryu (2000). The dynami- 
cal timescale associated with the remnant is given by 



1 / R \ V2 f 



M 



-1/2 



(60) 



the time it takes the binary remnant to achieve equilib- 
rium following coalescence. This is also the time that 
it takes the star, if it is driven far out-of-equilibrium 
by magnetic braking, to undergo collapse. The central 
rotation period of the remnant is 



tr 



R 



2ir ( 

sa 3 

ft ' V 20km 



3/2 



/ M 



-1/2 



(61) 



V3M Qy 

while the period at the equator is about three times 
longer. The timescale for magnetic braking of differential 
rotation by Alfven waves is given by 

-1/2 , M s 1/2 



Bn 



io 12 Gy V 2 0kmy V 3M © 



*A = — * 10 2 

VA 

(62) 

On this timescale, the angular velocity profile in the star 
is significantly altered. Finally, viscous dissipation drives 
the star to a new, uniformly rotating equilibrium state 
on a timescale 

^ 9 f R \ 23/4 f T \ 2 ( M ^ 5/4 

(63) 

where 77 = 347/j 13 / 4 r^ 2 cm 2 s" 1 (Cutler & Lindblom 
1987). In our simulations we choose parameters to pre- 
serve the inequality id yn < trot < tA < t v , although the 
relative magnitudes are altered for numerical tractability. 

In section l5.ll we will discuss results of dynamical sim- 
ulations with n — 0.001 and compare with the results of 



/ - 2 x 10 y ( — — 
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Paper I. Then in section 15.21 we treat n = 3 cylinders 
in detail as these case s di splay typical compressible be- 
haviors. Sections 15.31 and 15.41 will briefly review results 
of dynamical simulations for n = 1 and n = 5, respec- 
tively. Appendix |0 summarizes our numerical method 
and difference equations. 

5.1. Nearly Incompressible Cylinders 

As a check of our code, we first reproduce the analytic 
results of Paper I for incompressible cylinders by consid- 
ering a small value of the polytropic index: n = 0.001, 
which corresponds to T = 1001.0. This extremely stiff 
equation of state is a good approximation to an incom- 
pressible fluid. The polytropic index n was not simply 
taken to zero because n = corresponds to infinite sound 
speed, c s . The timestep in our simulations is limited by 
the Courant stability criterion St < Sr/cg, where Sr is 
the smallest spatial grid interval. Thus our numerical 
method is not applicable to the strictly incompressible 
n = case. 

Results for Case I (n = 0.001) with zero shear viscosity 
(?7 = 0.0) are shown in Figs. [21 and [21 In Fig- El we track 
the evolution of various energies with time and observe 
how the contributions to the total conserved energy os- 
cillate between rotational kinetic energy of the fluid and 
the energy of the azimuthal magnetic field. The angu- 
lar momentum is also strictly conserved. The period of 
these oscillations, Pa = 0.82 t& — 1.64 also agrees with 
the analytical result. We will henceforth define 



P A = 0.82 t A 



(64) 



so that the actual value of the timescale Pa will depend 
on the parameters of the model in question. The dynam- 
ical time in this case, idyn = 0.13, is much shorter than 
the Alfven wave oscillation period. Thus, the conver- 
sion of rotational to magnetic energy occurs over several 
dynamical timescales. We note that the outer radius 
of the cylinder remains constant to around one part in 
10 4 , showing that the effects of the slight degree of com- 
pressibility are not important for this simulation. There 
is essentially no radial motion and the cylinder remains 
very close to virial equilibrium. Let the normalized virial 
sum be defined as 



V = 



2U + 2T + W 

W\ ' 



(65) 



since \W\ is constant (see eqs. ^7j and ^Bl)- F° r 

this evolution, V < 1.5 x 10~ 5 (V = for equi- 
librium configurations). As an additional check, we 
found that the rotational and magnetic energies at 
times t = Pa/4, 3Pa/4, . . . when the azimuthal magnetic 
field reaches its maximum magnitude, are as follows: 
Pkin/Pki„(0) = 0.5136 and P mag /P ki n(0) = 0.4865. 
These agree well with the analytical results given in Pa- 
per I of 0.513 and 0.487, respectively. 

A crucial property of the analytic solution for the case 
n = is the scaling behavior (see Paper I) . For that solu- 
tion, the amplitudes of all evolved quantities are entirely 
independent of the initial radial seed field given in equa- 
tion (|34|l . Only the timescale of the oscillations of the 
energy components depends on B and this timescale is 
proportional to Bq . We confirmed that our numerical 
results for n = 0.001 also satisfy this scaling behavior 



h >\ 



0.8 



a 

H 0.4 -/ 



0.2 



Fig. 2. — Energy evolution for a differentially rotating, nearly 
incompressible star with zero viscosity and a vacuum exterior 
(Case I). The dashed line shows i?km and the dot-dash line shows 
Emag (E grav and Ej nt are essentially constant and are not shown). 
All energies are normalized to i?kin(0). The sum of -Ekin and E mag 
is conserved and remains equal to the initial rotational energy 
^kin(O), which is plotted as the solid horizontal line at E nerg y 
= 1. Time is in nondimensional units according to equation 1531 . 



to high accuracy. This check was performed by dou- 
bling Bq and observing that the oscillation period halves: 
Pa — * Pa/2. This scaling is not expected to hold for 
general values of n, however, because input parameters 
characterizing the initial rotation and magnetic field en- 
ter explicitly in the evolution system (see eqs. |55j-|58j). 
However, the analytic scaling should be approximately 
valid for stiff equations of state in a first approximation. 

The tradeoff between rotational kinetic and magnetic 
energies, Pkin and E mag , seen in Fig. [21 matches the 
standing Alfven wave behavior seen in Paper I. First the 
differential rotation generates a nonzero B^ which drains 
energy away from the differential rotation. When B^ 
reaches its maximum, the cylinder is uniformly rotating, 
and the built-up magnetic stress begins to drive differ- 
ential rotation in the opposite direction, unwinding the 
magnetic field and then winding it up again in the oppo- 
site sense. This cycle corresponds to a standing Alfven 
wave with period Pa- This process is shown explicitly 
in Fig. [3] which gives cross-sectional views of represen- 
tative magnetic field lines at critical phases during an 
oscillation period. The field lines were created by inte- 
grating the equation rd<p/dr — B^jB r to find <f> = 4>(r). 
(Our 1+1 Lagrangian evolution determines B^ and the 
Eulerian coordinate r for a given fluid element, while B r 
is fixed.) At times Pa/4 and 3Pa/4 corresponding to 
the maximum magnitude of B^, the field lines are highly 
twisted. The degree to which the field lines are twisted 
depends on the magnitude of the initial radial magnetic 
field (see eq. [2H] an d associated discussion). Even for a 
small initial field, however, the azimuthal magnetic field 
will grow to the same high value sufficient to brake the 
differential motion and drive the star to oscillate about 
the state of uniform rotation. This will have important 
consequences for a hypermassive star which depends on 
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Fig. 3. — Magnetic field line configurations for Case I (nearly 
incompressible, n = 0.001, zero viscosity) at selected times. The 
twisting of the field lines is due to differential rotation. Snapshots 
are given at time inte rval s of Pa/4, where Pa is the standing Alfven 
wave period (see eq. I64l t. Each frozen-in field line passes through 
the same fluid elements at all times. The bold field line is an 
arbitrarily chosen fiducial line. 



Fig. 4. — Energy evolution for Case I (n = 0.001) with shear vis- 
cosity given by r\ = 0.2. The dotted line labeled "iSgrav" is actually 
(-Egrav — S grav (0))/Pkin(0). The line labeled "Ei nt " is also defined 
in this way. The expansion of the cylinder causes the gravitational 
potential energy to become less negative, compensating for the loss 
of -Ejjin ■ The sum of the energies is conserved and is plotted as the 
horizontal solid line labeled "E total ." 



differential rotation for stability against gravitational col- 
lapse. 

Next, we consider an approximately incompressible 
model with nonzero shear viscosity (Case I with 77 = 0.2), 
and results are given in Figs. 0] and [3] In Fig.^J we track 
the evolution of the various energies with time. We ob- 
serve how the oscillations of the rotational kinetic and 
azimuthal magnetic field energies are now damped by 
viscosity. The damping takes place over a few periods 
(Pa) of the magnetic field oscillations, consistent with 
the choice of 77 = 0.2 and hence t v /tA — 2.5. This 
damping behavior was also seen in the results of Pa- 
per I. However, in this case, the heating due to viscosity 
causes the cylinder to expand, with the outer radius ex- 
panding by about 7%. The expansion causes a decrease 
in the magnitude of the (negative) gravitational poten- 
tial energy .Egrav which compensates for the loss of rota- 
tional kinetic energy due to the viscous dissipation. (We 
note that the line labeled "Egrav" in Fig. 01 is actually 
(Eg r av — E g rav(0))/Ekin(0).) In contrast, the internal en- 
ergy does not change much from its initial value. We 
note that total energy and angular momentum are con- 
served very well for this simulation (to one part in 10 7 
and 10 9 , respectively). The fact that the cylinder ex- 
pands in this case clearly indicates that, when there is a 
significant amount of viscous heating, compressibility ef- 
fects are significant even for n = 0.001. To confirm that 
this expansion is indeed caused by heating, we performed 
this simulation a second time without the viscous heating 
terms (i.e. by removing the terms proportional to 77 in eq. 
|58j1. No expansion occurs in this case. When viscous 
heating is included, the pressure term II in the virial equi- 
librium equation (eq. |17|) significantly increases. This 
change is not reflected in the internal energy, however, 
since Ei nt = nil (by eqs. ^B] and [SD| and the T-law 
EOS) and n is small in this case. 



5.2. The n = 3 Results 

We now present results for the slowly differentially ro- 
tating n — 3 case (Case II) without viscosity, and with 
(3 = 0.030 (< /3 max ). The numerical evolution for this 
model shows a slight overall expansion of the cylinder and 
radial oscillations which increase in magnitude toward 
the surface. Oscillations with the same period are also 
seen in the radial profile of B,p, as shown in Fig. [HI which 
shows snapshots of the profile at selected phases of 
the quasi-oscillation period. This period is ~ 0.244 Pa, 
significantly shorter than the Alfven period, Pa, as de- 
fined in equation Q64[l. The dynamical timescale for this 
case, tdyn = 0.0091 is much shorter than the period of 
the radial oscillations (since 0.244 Pa = 0.400), and the 
oscillatory behavior is therefore quasistatic. The angular 
velocity profiles for this model also display this oscilla- 
tory behavior, with the angular velocity becoming more 
uniform at phases corresponding to maxima in P^. The 
radial oscillations show no indication of decaying during 
the simulation, which ran for several Alfven timescales: 
tfinal ~ 3£a = 3.7 Pa- Thus, with no dissipation mecha- 
nism, these oscillations presumably continue indefinitely. 
[Note, however, that a small amount of viscosity was 
added in order to stabilize the run, so that = 70.0 1&. 
This has no significant impact on this run, since the total 
length of the run was only ~ 3 *a-] That this model does 
not display any strong radial motions is expected since 
= T/\W\ is below the limit for uniform rotation and 
hence the configuration can be driven to uniform rota- 
tion quasi-periodically without undergoing appreciable 
contraction and/or expansion. 

Now suppose a significant shear viscosity is introduced 
to the case described in the previous paragraph. As in 
Case II without viscosity, the cylinder expands slightly. 
Radial oscillations occur again with roughly the same pe- 
riod (P w 0.244 P4), but these oscillations are damped 
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Fig. 5. — Azimuthal magnetic field profiles for Case I with viscos- 
ity. The radius r and azimuthal field are plotted in nondimen- 
sional units according to equation 1531 . Snapshots of the profiles 
are taken at the same times as in Fig. [3] The behavior is ap- 
proximately periodic, with a period similar to Pa- Because of the 
viscosity, oscillations of the profile are damped. Note that the 
profiles do not meet at the radius of the outer mass shell due to 
the slight expansion of the cylinder. 
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Fig. 6. — Azimuthal magnetic field profiles for Case II (n = 3 
with slow differential rotation) for r\ = 0. The radius r and az- 
imuthal field B$ are plotted in nondimensional units according 
to equation 1531 . Snapshots of the profiles are taken at steps of 
0.061 Pa- The behavior is seen to be quasi-periodic, with a pe- 
riod of ~ 0.244 Pa- For example, the maximum magnitudes of B^ 
occur approximately at fractions 1/4,3/4, and 5/4 of the period. 
This is the same period seen in the gentle radial oscillations which 
occur during the evolution of this model. 



after a few periods due to the shear viscosity. The final 
state is uniformly rotating with no azimuthal magnetic 
field. In this case, the cylinder easily relaxes to uniform 
rotation and dynamical equilibrium within a few oscilla- 
tion periods of the magnetic field. The excess rotational 
energy is converted into internal energy (heat). 

We now examine the results for the rapidly differen- 
tially rotating n — 3 case (Case III) with -q = and 
(3 = 0.191 3> Anax- Though we are interested in the 
behavior with no shear viscosity, we again add a very 
small shear viscosity to stabilize the simulation. In this 
case, taking r\ such that — 500.0 1& is sufficient, and 
the effects of the shear viscosity will not be seen in the 
following results. In all cases discussed hereafter, we will 
add a similar small, stabilizing viscosity that does not af- 
fect the results on the timescales of interest. (This small 
viscosity will, of course, not be added in cases where we 
are considering the effects of a significant shear viscos- 
ity with 77 chosen such that t n — 2.5 The motion 
of the Lagrangian mass tracers for this case is shown in 
Fig. [7f a) (Fig.[7{b) shows the analogous results with sig- 
nificant shear viscosity and will be discussed below). The 
inner ~ 95% of the mass undergoes a significant contrac- 
tion and a slight bounce. The inset figure shows that the 
outer ~ 5% of the mass expands to large radii, forming 
a diffuse atmosphere. In the final state, the interior of 
the star is slowly rotating, having lost most of its angular 
momentum to the outer envelope. 

Figure 12a) shows that the duration of the contrac- 
tion is approximately the Alfven timescale, £a = 2. 
The dynamical timescale as defined in equation (|60[) is 
much shorter than the Alfven timescale: £a = 66.4 t^ yn . 
Since the contraction takes place over many dynamical 
timescales, the cylinder is never far from equilibrium and 
the evolution is quasistatic. The evolution of the virial 



sum (eq. [123); which is shown in Fig. [8] illustrates this 
fact. While the pressure and kinetic energy terms of the 
virial sum depart significantly from their initial values, 
the sum V remains near its equilibrium value of zero. 
Also, typical radial velocities of the contracting shells 
have magnitudes only ~ 10% of the "free-fall velocity," 
Uff, where vg = R/td yn - I n & realistic relativistic hy- 
permassive star, this quasistatic contraction would likely 
lead to dynamical radial instability and then catastrophic 
gravitational collapse of the core. The diffuse atmosphere 
seen in our simulation suggests that magnetic braking op- 
erating in rotating stars may also lead to the ejection of 
winds or a diffuse ambient disk. 

The time evolution of the various components of the 
energy is shown in Fig.El The most salient feature is the 
strong contraction, which results in the sharp increase 
in £"i n t and an increase in the magnitude of the gravi- 
tational potential energy. One also sees that, early in 
the simulation, the magnetic energy grows to an appre- 
ciable fraction of the initial rotational kinetic energy. In 
the final state, most of the kinetic energy has been con- 
verted into internal energy by the contraction. Angular 
momentum is nonetheless conserved because the outer 
shells are very slowly rotating at very large radii. We 
show cross-sectional views of the magnetic field lines in 
Fig. 1101 This shows that the maximum magnitude of 
occurs early in the run, with some small oscillations fol- 
lowing. There is a striking contrast in behavior between 
this case with rapid differential rotation and Case II dis- 
cussed above. For Case II, (3 < /3 max and the magnetic 
field (and viscosity, when present) can drive the cylin- 
der toward uniform rotation without significantly alter- 
ing the structure of the cylinder. In Case III, however, (3 
is significantly larger than /3 max and the cylinder cannot 
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Fig. 7. — (a) Spacetime diagram featuring Lagrangian mass tracers for Case III (n = 3, rapid differential rotation) with r\ = 0. The 
solid lines show mass fractions 0.10, 0.20, . . . , 1.0 while the dotted lines enclose fractions 0.91, 0.92, . . . 0.99. The large plot demonstrates 
the contraction of the inner shells and the expansion of the outer envelope. The inset shows the same outermost tracers out to a much 
larger radius and for a longer time. This inset shows that a diffuse atmosphere has been formed outside of the cylinder, (b) Same as (a) 
except now with viscosity r) = 0.2. The contraction of the inner shells and expansion of the outer envelope shown here are milder than the 
corresponding behaviors in (a) with zero viscosity. 
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Fig. 8. — Individual components of the virial sum V = (2II+2T+ 
W)/\W\ioi Case III (V = in strict equilibrium) with r) = 0. The 
components 211, 2T, and W are normalized to \W\ and are plotted 
with their initial values subtracted off. That V remains small while 
the cylinder undergoes significant contraction indicates that the 
contraction and readjustment occur in a quasistatic fashion. 
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Fig. 9. — Energy evolution for Case III (n = 3 and rapid differ- 
ential rotation) with r\ = 0. All energies are defined as in Fig. [I] 
The rapid rise in E; nt and sharp drop in E gISjV correspond to the 
contraction of the core. The magnetic energy -E m ag grows to an 
appreciable fraction of the initial kinetic energy before decays 
back to zero in the final state. 



relax to uniform rotation without restructuring. In this 
case, we see both a strong contraction and expansion of 
the outer layers. Thus, when the differential rotation is 
strong, magnetic braking leads to dramatic changes in 
the stellar structure. 

These results can be compared with the analogous case 
with significant shear viscosity present (ry = 0.2). Fig- 
ure [7f b) shows the Lagrangian matter tracers for this 
case. A contraction similar to that seen in Case III with 



r\ = is seen again here. However, the inset plot shows 
that the outer shells are not ejected to such large dis- 
tances as they are in the previous case. Also, a strong 
bounce is not seen in the present case. One may conclude 
that the presence of viscosity has resulted in slightly 
milder behavior. This is also seen in the fact that the 
central density increases by a factor of 3.50 in this case, 
but 3.71 in the case without significant shear viscosity. 
The approach to uniform rotation in this case is easily 
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Fig. 10. — Magnetic field line configurations for Case III (n = 3, 
rapid differential rotation) with r) = as in Fig. [3] The magnetic 
field configuration at each time is shown only for the inner 96% 
of the mass because the outer envelope is ejected to large radius, 
forming a low-density atmosphere (hence the continuation of each 
field line outside the cylinder is not shown). The azimuthal mag- 
netic field is clearly strongest at early times. 
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Fig. 11. — Angular velocity profiles for Case III with r) = 0.2 at 
selected times. The configuration is driven toward slow, uniform 
rotation. Later profiles extend to r > 1 because of the expan- 
sion of the cylinder. Angular velocity and tim e are expressed in 
nondimensional units according to equation 1531 . Though the final 
uniform value of Q is much smaller than the initial values in the 
interior, angular momentum is conserved due to the slow rotation 
of the outer shells at large radii. 



seen in Fig. 1111 which gives angular velocity profiles at 
various times. As before, the final state has very little 
rotational kinetic energy and angular momentum is con- 
served by transferring angular momentum (via viscosity 
and magnetic fields in this case) to slowly rotating outer 
shells at large radii. 

The final state for Case III with r\ = has 
(P/p r )/(P/p r ) t= a = 1.11, which indicates that the 
strong contraction and bounce resulted in significant 



shock heating. [The symbol (P/p T ) represents an av- 
erage of P/p T over the inner 96% of the mass. The 
average was restricted to the inner 96% because P/p r 
could not be calculated accurately in the low-density at- 
mosphere.] The existence of shocks is not in conflict 
with the overall quasistatic nature of the evolution be- 
cause all of the shocks that we have observed are weak 
(with Mach numbers in the range 1.0 < M. < 2) and 
occur in the outer, low density regions. Consider the 
shock which formed during this evolution at t ~ -Pa/4 
at p(r)/p t ~ 0.93. The MHD shock propagated outward 
into the atmosphere of the cylinder, and the incoming 
fluid had a slightly supersonic inward radial velocity with 
respect to the shock front. The angular velocity and az- 
imuthal magnetic field had discontinuously larger values 
at the front in the unshocked fluid versus the shocked 
fluid. This is because the shells falling into the shock 
front spin up due to momentum conservation. The re- 
sulting discontinuity in f2 then leads to a discontinuity 
in B^ due to flux freezing. (In Section 15.41 we will give 
an example of a shock for which actually changes 
sign across the shock front.) We explored the extent to 
which the quantities upwind and downwind of this shock 
satisfy the jump conditions for shocks in magnetic fluids 
(Landau, Lifshitz, & Pitaevskii 1984). For example, the 
following condition can be derived from continuity of the 
z— component of the electric field (which applies for a 
steady state shock) and conservation of mass across the 
MHD shock front: 

/ -B</>,2 -60,1 \ B r 
Pi = {v<f,,2-v</, t i), (66 

V P2 Pi J V r ,l 

where the subscripts "1" and "2" refer to the pre- and 
post-shock fluids respectively (this equation is in physical 
units, not nondimensional units). The difference between 
the right and left hand sides of equation (|66|l for the 
shock in our numerical solution is ~ 8%. The other jump 
conditions are also satisfied with errors of < 20%. This 
agreement is reasonable considering the simplicity of the 
artificial viscosity scheme employed by our code and the 
fact that the jump conditions with which we compared 
our results assume steady state rather than dynamical 
conditions. 

5.3. The n = 1 Results 

We now consider the differentially rotating n = 1 case 
(Case IV) without viscosity. This case has (3 = 0.168, 
which is below the upper limit for uniform rotation, /3 ma x- 
Recall that all differentially rotating n = 1 models built 
according to equation ifFf)! have (3 < (3 max . However, 
(3 = 0.168 is near the upper limit of (3 for the given ro- 
tation law. Even though this case can relax to uniform 
rotation without appreciable restructuring, we find that 
this does not happen. The star contracts slightly so that 
the central density grows by a factor of 1.819 and sets up 
a large diffuse atmosphere (see Fig. 11211 . FigurelT3lshows 
the evolution of the different contributions to the energy. 
We find that the energy in magnetic fields increases to a 
significant fraction of the initial kinetic energy, but most 
of the initial kinetic energy is ultimately converted into 
internal energy through shock heating. It appears that 
the final state will be a hot, uniformly rotating star with 
a dense core and a diffuse atmosphere. We can also see 
oscillations on the energy plot with a period of ~ 0.76 Pa- 



14 




0.5 1 1.5 



r 

Fig. 12. — Density profiles for Case IV (n = 1) with zero viscosity 
(r; = 0). The radius r and de nsity p are plotted in nondimensional 
units according to equation 1531 . This plot shows a contraction 
of the inner material with the outer material forming a diffuse 
atmosphere. 



This is close to the Alfven period, Pa, as defined in equa- 
tion (|64|) . The dynamical timescale for this evolution is 
£dyn = 0.111, which is much shorter than the Alfven 
time. The virial diagnostic shows that this evolution is 
quasistatic and similar to the n = 3 cases. We now con- 
sider the effects of adding a shear viscosity. As before, 
the cylinder develops a large diffuse atmosphere and con- 
tracts, but does so slightly less severely than in the case 
with r\ = 0. Radial oscillations occur again with roughly 
the same period (P f=s 0.76 Pa), but these oscillations are 
damped after a few periods due to the shear viscosity. In 
the n = 1 inviscid case, the oscillations make determin- 
ing the final state difficult, but, as can be seen in Fig. 
1141 the cylindrical model with viscosity tends to uniform 
rotation by a time t ~ Pa- 

5.4. The n = 5 Results 

We now describe results for the slowly differentially ro- 
tating n — 5 case (Case V) without viscosity. For this 
case, we choose (3 = 0.0183 < f3 max . Therefore, we do 
not expect strong radial motions in this case. This case 
behaves much like the analogous n — 3 case (Case II), 
so we do not show any plots for it. The major difference 
is that this cylinder experienced a slight contraction and 
radial oscillations which increase in magnitude toward 
the surface. The central density grows only by a factor 
of 1.170 over three Alfven periods. We also see oscil- 
lations in the magnetic field with a period ~ 0.195 Pa, 
much shorter than the Alfven period. The dynamical 
timescale for this evolution is tdyn = 0.0020, which is 
again much shorter than the Alfven time. As in the 
n — 3 cases, the virial diagnostic shows that this evo- 
lution is quasistatic. The only mechanism to dissipate 
energy is through shock heating, and, since the behavior 
is mild, this would require a very long run to reach a 
final state. Adding a significant shear viscosity to this 
model results in damping of the radial oscillations and 
rapid approach to uniform rotation, as in Case II. 
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Fig. 13. — Energy evolution for Case IV (n = 1) with r) = 0. 
All energies are defined as in Fig. [3] The sharp increase in Si n t 
and sharp decrease in E gIav correspond to the contraction of the 
interior shells. Then as the outer layers continue to expand, -Egrav 
increases. Also, -E m ag grows to an appreciable fraction of the initial 
kinetic energy before approaching zero. Note the oscillations in the 
energy components. The sum of the energies is conserved and is 
plotted as the horizontal solid line. 
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Fig. 14. — Angular velocity profiles for Case IV (n = 1) with 
7] = 0.2. The radius r and angular vel ocity Q are plotted in nondi- 
mensional units according to equation 1531 . Viscosity and magnetic 
braking drive the star to uniform rotation. 



Next, we examine the results for the rapidly differen- 
tially rotating n = 5 case (Case VI) with 77 = and 
where (3 = 0.361 ^> /3 max - The motion of the Lagrangian 
mass tracers for this case is shown in Fig. The inner 
~ 91% of the mass undergoes a severe contraction. The 
inset figure shows that the outer ~ 9% of the mass ex- 
pands to large radii, forming a diffuse atmosphere. The 
time evolution of the various components of the energy is 
shown in Fig. 1161 which clearly shows the strong contrac- 
tion. As in the n = 3 cases, we see that magnetic braking 
leads to significant changes in the structure when the ro- 



MAGNETIC BRAKING OF DIFFERENTIAL ROTATION 



15 



5 



t. 




12 3 

r 



Fig. 15. — Spacctimc diagram featuring Lagrangian mass tracers 
for Case VI (n = 5) without viscosity (rj = 0.0). In the large 
plot, the solid lines enclose mass fractions 0.10, 0.20, . . . , 1.0 while 
the dotted lines enclose fractions 0.92, 0.94, 0.96, and 0.98. In the 
inset, the solid lines have the same meaning, while the dotted lines 
enclose fractions 0.91,0.91 0.99. The outermost line in the 
inset encloses a mass fraction of 1.0. The growth of the azimuthal 
magnetic field causes a severe contraction of the inner shells and 
large ejection of the outer layers. The inset shows that the outer 
layers form a diffuse atmosphere. 
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Fig. 16. — Energy evolution for Case VI (n = 5) with r) = 0. All 
energies arc defined as in Fig. |3] The sharp increase in Si n t and 
sharp decrease in E gTav again correspond to the contraction of the 
interior shells. Also, E ma g grows to an appreciable fraction of the 
initial kinetic energy before approaching zero. 



tation is strong. The final state for Case VI with 77 = 
has (P/p r ) — 2.65, which indicates that the strong con- 
traction resulted in significant shock heating as in Case 
III. FigureElgives a snapshot of an outgoing shock prop- 
agating through the envelope in the Case IV evolution. 
The grayscale shows [P/ p r )/[P/ p r ] t= o, a local measure 
of heating. The solid lines represent magnetic field lines, 
which are clearly discontinuous at the shock front. How- 
ever, examining the virial again reveals that this evolu- 
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Fig. 17. — Outgoing shock for Case VI with rj = at time 
t = 0.24 Pa = 0.40. (Only one quadrant is shown.) The shock front 
encloses ~ 83% of the total mass. The ratio [P/p r ]/[P/p r ]t— is 
given by the gray shading, with magnitudes shown in the bar on 
the upper right. The solid lines are representative magnetic field 
lines, calculated as in Fig. [3] The discontinuity indicates a change 
in the sign of B^. 
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Fig. 18. — Spacetime diagram featuring Lagrangian mass tracers 
for Case VI (n = 5) with v isco sity (77 = 0.2). The lines have 
the same meaning as in Fig. 1151 Notice that the contraction is 
slightly less severe and the ejection much less severe than the case 
without viscosity. 



tion is nearly quasistatic. These results can be compared 
with the analogous case with significant shear viscosity 
present (77 = 0.2). Figure ITS1 shows the mass tracers for 
this case. A contraction similar to that seen in Case VI 
with rj — is seen again here. However, the inset plot 
shows that the outer shells are not ejected to such large 
distances as they are in the previous case. Also, the cen- 
tral density increases by a factor of 135.9 in this case, 
but 153.39 in the case without significant shear viscos- 
ity. One may conclude that the presence of viscosity has 
again resulted in slightly milder behavior. 
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Table 3 
Summary of Results 



Initial Data Characteristics 


TlA/n Pi TYi iiTi 1 Ron Q\ri(~if 

L_s y lld.lllH_.Ojl 1—) <_ Lid' V 1U1 


Final State Characteristics 8, 


Case 


n 

p 


n 


V 




(P/p r ) (R/Roho 




Pc/pc(0) 


I 


0.043 


0.001 


0.0 


crt" a ti H i n cr Alfvpn wavp 

r3 LiCLllU-lllt ill! V (jll VVOiVIj 


1.00 


1.00 


1.00 


1.00 








0.2 


U.CLlllJJCLI. Uol^llld L 1(JI1D 




1.02 


1.08 


1.00 


II 


0.030 


3.0 


0.0 


PYtifln t;i r>n riQPi 1 1 n t i r>n q 

CA.JJcL11oHJ11 , USL-lllCtLlUllo 


1.00 


0.99 


1.00 


1.02 








0.2 


slight expansion 


1.00 


0.99 


1.00 


1.02 


III 


0.191 


3.0 


0.0 


core contraction 


1.11 


0.52 


0.57 


3.71 








0.2 


core contraction 


1.40 


0.53 


0.63 


3.50 


IV 


0.168 


1.0 


0.0 


slight contraction 


1.07 


0.80 


0.88 


1.82 








0.2 


slight contraction 


1.22 


0.82 


0.95 


1.77 


V 


0.018 


5.0 


0.0 


slight contraction, oscillations 


1.00 


0.92 


0.91 


1.17 








0.2 


slight contraction 


1.00 


0.92 


0.91 


1.17 


VI 


0.361 


5.0 


0.0 


strong core contraction 


2.65 


0.08 


0.78 


153.4 








0.2 


strong core contraction 


2.54 


0.09 


0.73 


135.9 



a (P/p v ) in units of (P/p r )t=o', values larger than 1 indicate heating. (R/Ro)i denotes the radius containing a fraction i 
of the mass, in units of its initial value. 



6. CONCLUSIONS 

Poloidal seed magnetic fields and viscosity have impor- 
tant dynamical effects on differentially rotating neutron 
stars. Even for a small seed magnetic field, differential ro- 
tation generates toroidal Alfven waves which amplify and 
drive the star toward uniform rotation. This magnetic 
braking process can remove a significant amount of rota- 
tional energy from the star and store it in the azimuthal 
magnetic field. Though it acts on a longer timescale, 
shear viscosity also drives the star toward uniform rota- 
tion. For a hypermassive star supported against collapse 
by differential rotation, magnetic braking and viscosity 
can lead ultimately to catastrophic collapse. 

The strength of the differential rotation, the degree of 
compressibility, and the amount of shear viscosity all af- 
fect the response of differentially rotating cylinders to 
the initial magnetic field. Our calculations have shown 
very different behavior when (3 — T/\W\ is below the 
upper limit for uniform rotation, /3 ma x, than when it 
is above this limit. Simulations for n — 3 and n = 5 
with j3 < /3max show that the cylinders oscillate and ei- 
ther expand or contract slightly to accommodate the ef- 
fects of magnetic braking and viscous damping. Large 
changes in the structure are not seen for these cases. 
For > /3 max , on the other hand, the outer layers are 
ejected to large radii, while most of the star contracts 
quasistatically. In a simulation of a relativistic hyper- 
massive star with more realistic geometry, this behavior 
would likely correspond to quasistatic contraction lead- 
ing to catastrophic collapse and escape of some ejected 
material. The dynamical behavior is more extreme for 
models with greater compressibility. This result is rea- 
sonable since softer equations of state allow stronger ra- 
dial motions. Inclusion of a shear viscosity often moder- 
ates the behavior of an evolving model, pacing angular 
momentum transport and damping the toroidal Alfven 
waves which arise due to differential rotation. In partic- 
ular, the numerical simulations in the rapidly rotating 
n = 3 and n — 5 cases show that the contraction of 
the core and ejection of the outer shells is milder when 
a significant shear viscosity is present. Our results are 
summarized in Table 3. 

Though our model for differentially rotating neutron 
stars is highly idealized, it accommodates magnetic 



fields, differential rotation, viscosity, and shocks in a 
simple, one-dimensional Lagrangian MHD scheme. In 
addition, we are able to handle the wide disparity be- 
tween the dynamical timescale and the Alfven and vis- 
cous timescales. This disparity will likely prove taxing for 
relativistic MHD evolution codes since it will be neces- 
sary to evolve the system for many dynamical timescales 
in order to see the effects of the magnetic fields. Be- 
cause the cylindrical geometry of the present calculation 
greatly simplifies the computational problem, this class 
of numerical models is useful in developing intuition for 
the physical effects which almost certainly play an impor- 
tant role in nascent neutron stars. More realistic evolu- 
tionary calculations of magnetic braking in neutron stars 
will serve to further crystallize our understanding of the 
evolution of differentially rotating neutron stars, includ- 
ing hypermassive stars that may arise through binary 
merger or core collapse, and may shed light on the phys- 
ical origins of gravitational wave sources and gamma ray 
bursts. Several additional physical effects are expected 
to be important. For example, Shapiro (Paper I) showed 
that the presence of an atmosphere outside of the star 
allows efficient angular momentum loss by partial trans- 
mission of Alfven waves at the surface. This will have 
important effects on the dynamical behavior of a neutron 
star relaxing to uniform rotation. In addition, our treat- 
ment does not account for evolution of the magnetic field 
due to turbulence and convection. However, even with- 
out these additional physical processes, our calculations 
reveal qualitatively important features of the effects of 
magnetic braking on the stability of neutron stars. We 
hope to continue to pursue these questions through more 
realistic computational investigations, including the ef- 
fects of general relativity. 
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APPENDIX 

THE VIRIAL THEOREM 

In this appendix, we closely follow the derivation of the virial theorem in Shapiro & Teukolsky (1983) Section 7.1. 
We will discuss the inviscid case, because we do not expect viscosity to play a dominant role in equilibrium stars. We 
start by dotting the position vector x into equation l|31fl . We then multiply through by p, integrate over a cylindrical 
volume of length L and radius R, and divide by L (we will work with quantities per unit length as we did in section 
3.3). Define T, H, and W as in section 3.3 and the moment of inertia per unit length as 



/= / pr'dA. (Al) 
Jo 

The left hand side of equation (|3"T)l then becomes 

lw- 2T - < A2 » 

The gravitational term is W. We will calculate the pressure term in detail to display a feature of this geometry. The 
pressure term is 

- i J x- VPd 3 x = -i J V • {Px)d 3 x+j J PV • (x)d 3 x 

= -n + 3n, (A3) 

where we have used 

=-n. (A4) 

The difference between this infinite cylinder case and the case for bounded distributions is that we must keep a term 
from the endcaps when applying the divergence theorem in equation (|A4I) . Applying similar manipulations to the 
magnetic terms shows that these terms cancel. Collecting these results gives 

1§ = 2T + W + 2U. (A5) 
Assuming steady-state (d/dt — 0) gives the virial theorem quoted in equation (|17|l . 

STABILITY OF CYLINDRICAL POLYTROPES 

In this section, we derive by means of an energy variational calculation the critical polytropic parameter r cr it for 
the stability of a slowly rotating cylinder against radial perturbations, adapting the treatment of stability for ordinary 
stars given in Shapiro and Teukolsky (1983). Because we are considering equilibrium models, there is no azimuthal 
magnetic field, and we can take the total energy as 

E = Pint + T + -Egrav , (Bl) 

where T is the rotational kinetic energy. Consider a sequence of rotating equilibrium cylinders with constant angular 
momentum per unit length, J, parameterized by central density, p c . The equilibrium mass per unit length p t is 
determined by the condition dE/dp c = 0. Stability then requires d 2 E/dp 2 > 0. Using the T-law EOS for polytropes, 

Pint = / d^jr^Y = kiKp t p r c ~\ (B2) 

where k\ is some nondimensional constant of order unity that depends on the density profile through Y for slowly 
rotating configurations. The kinetic energy is just T = J 2 /2I, where I is the moment of inertia per unit length. Since 
/ cx p,tR 2 oc (Lt|/p c , one may write 

T = k 2 ^, (B3) 
K 

where k 2 is another constant. (We note that ki and k 2 may be derived analytically for nonrotating cylinders by 
rewriting the integral in terms of Lane-Emden variables.) 
Next consider the gravitational potential energy. As in equation (|50|) , the gravitational energy per unit length is 

£ g rav = \ J dr2irrp<S> = \J dr^$. (B4) 

Now from equation JSJl, d<&/dr = 2p(r)/r. Integrating this equation inward from a fiducial radius r outside of the 
cylinder gives 

$(r) = 2p t ln(r/r ), (r > R). (B5) 
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We now integrate the right-hand side of equation l|B4|) by parts, insert the forms for <f> and its radial derivative, 
and obtain ^ 

£ grav = HR/ro) - I dr^±. (B6) 



Any terms indepen dent of p c may be discarded for the purposes of this variational calculation. The second term on the 
r. h. s. of equation l|B6|) scales as p 2 and is independent of p c . However, the first term depends on p c since R 2 oc pt/p c - 
Then 

-Egrav = \p% hi I j ) + constants. (B7) 

From equations (|B2|) . (|B3|) . and (|B7|I . the total energy is (up to additive constants) 

E = k 1 p t Kp r c - 1 + k 2 ^ + ^ t \n(-^ 1 ). (B8) 

Mi \Pc r o/ 

The conditions for equilibrium and for the onset of instability become 



ap t2 ,.2 



2 =o=fc 1 (r-i)(r-2) Mt A>r j + f^- (B9) 

The first equation is equivalent to the condition for hydrostatic equilibrium (see eq. |10|). Inserting the first equation 
into the second to eliminate the term proportional to fci gives, after rearranging, 

i (r -i)-(r- 2)^-i=o. (bio) 
Pt Pt 

Now, by equa tion ljl8|l . Pt = \W\, where W is the gravitation al ter m appearing in the virial theorem. We also use 
equation (|B3|> to identify the kinetic energy T. Since equation (|B10|) is the condition marking the onset of instability, 
solving for T gives the critical value: 

r °* = [! = Z/m (cylinders) - (B11) 

By contrast, the result for ordinary bounded stars is (Shapiro & Teukolsky 1983) 

r CT it = - V ~ 5T ' ,n^) (bounded stars). (B12) 
3 (1 - 2T/\W\) v ' ' 

Note that, for nonrotating cylinders with T/\ W\ = 0, the critical polytropic parameter is r cr j t = 1. For larger values of 
T/|W<^|, r cr j t decreases. But since T = 1 is the smallest T obtainable for any positive polytropic index n, all equilibrium 
cylindrical polytropes are radially stable. 

NUMERICAL TECHNIQUES 

We solve the MHD evolution equations using a one-dimensional Lagrangian finite differencing scheme in which the 
cylinder is partitioned in the radial direction into N shells of equal mass per unit length. The evolution equations are 
differenced on staggered meshes in time and space. The quantities r, v r , j, and p are defined on the shell boundaries 
while P,V= 1/p, e, B<p, and q, the artificial viscosity, are defined at the shell centers. Integral and half-integral spatial 
indices correspond to shell boundaries and shell centers, respectively. For example, one has where i = 1, 2, . . . , N + 1 
and -Pi+1/2 where £ = 1,2,..., JV. Updating at each timestep occurs in leap-frog fashion, with v r , q, and j defined at 

times with half-integral indices (< n+1 / 2 ) and r, P, V, B^, and e defined at times t n . 

Shocks are handled using artificial viscosity as described by Richtmyer and Morton (1967). In physical units, the 
artificial viscosity prescription is 

^m 2 , Xdv r /dr<0; (C1) 
0, otherwise. 

where 5r is the typical spatial mesh size and the thickness of the resulting smoothed shock is approximately a 2 5r 
(we usually take a = 1.75). The timestep is governed by the usual Courant condition in the absence of viscosity: 
St = b5r/c s , where c s is the local sound speed and b is an constant less than unity (we usually take b — 0.3). If shear 
viscosity is present, the condition becomes 

<5t = min(^,^l, (C2) 



c, 2D, 



off 



where D e g = 4rj/3p is an effective diffusion coefficient. 
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We found it convenient to use a different nondimensionalization scheme for our numerical scheme than was employed 
in the discussion above (eq. [3]). We will now use the following definitions: 



f = r/R, t = t/(R/c s (0)), fj = r,/(p c c a (0)R), (C3) 

j = j/(c s (0)R), B^B^/Bo, 

P = p/p c , v r = v r /c s (0), e = e/[Kp T c - 1 /{T-l)\ 

~p = p/( Pc R 2 ), P = P/(K P r c ), 

where c s (0) = (TKp^ -1 ) 1 ^ 2 is the sound speed at r = and t = 0. This differs from the previous scheme in that it relies 
on the sound speed a s a v elocity scale rather than the Alfven speed. We finite-difference the evolution equations 1)39(1- 
<|42l) , using equations (|C3|) to write the system in nondimensional form. (Below we drop the tildes on all nondimensional 
quantities.) If no shear viscosity is present, all derivatives are approximated by second order centered differences, so 
that the scheme is second order in space and time (except at the boundaries). The terms which account for shear 
viscosity are not time centered (i.e. these terms are first order in time). This simplification does not significantly affect 
our results since the viscous terms are small compared to the other terms in the evolution equations. 

Before giving the evolution system, we will introduce several definitions and conventions to simplify the notation. 
First, we will drop the r index on v r , the tj> index on and the tildes. We also define 8 = va/c s (0). In all of the 
following equations, the operator A, which takes a spatial difference, will be defined by AQi = Qi — Qi-i- In some of 
the equations, it is useful to define space and time averages for grid positions with half-integral indices. For example, 



i+l/2 



(r« +1 + r?)/2 and rf +1/2 
(eq. 001) becomes 

.ri+l/2 .n-1/2 



5t 



2tt^ 
Sp 



A(r; 



r™)/2. Then the evolution equation for specific angular momentum 
(r? +1 ) 3 + (rF) 3 ^ I '^ 1/2 



+ l/2-°i+l/2 



A? 



i+l 



A 



(C+i) 2 



(C4) 



We have used the fact that the Lagrangian mass increment between two shells, dp is a constant. The radial velocity 
is updated next as follows: 



v 



n+l/2 



n-1/2 



St 



2nrl 
TSp 
1 



A(iT +1/2 ) - 

+ 1/2^2 



2vr 



nO 2 
5pr\ 



-A((r7 +1/2 ) 2 5IV 1/2 ) 



2(C) 3 

87T77 

' 36 fir? ' 



A 



(c+i) 3 



Ar™ 



J i+1 



-1/2^ 



'i+l 



(C5) 



The appearing in the gravitational potential term of the above equation is just the nondimensional surface radius 
from the solution of the Lane-Emden equation. It arises here only because of our choice of nondimensionalization. 
Next, we update the locations of shell boundaries according to 

„"+l „n 



n+l/2 



St 



The specific volume of each shell is then updated to reflect the new shell boundaries: 



yn+l _ yO 



12 \ A(r? +1 ) 2 



(C6) 



(C7) 



where the superscript "0" refers to the initial value (this equation merely expresses the fact that the volume of a thin 
shell is dV = irLdr 2 ). We next integrate the evolution equation for the magnetic field. 



B 



i+l/2 



B n , 



1/2 



St 



i+l/2 



r>n 
B i+1 



Av 



n+l/2 - 



1 



/2) 



n+l/2 



'i+l 



n+l/2 
r r+l 



A 



.n+l/2 s 
3 i+l 

■ "+1/2^2 



(C8) 



This equation can be solved for the updated magnetic field in each cell, BV^, 2 . The next quantity to update is the 
artificial viscosity: 



n+l/2 
1i+l/2 



2g2r (Av n+1/2 ) 2 

+ 1/2 V L * U i+l ) > 



o, 



if Av^<0; 
otherwise. 



(C9) 



The specific internal energy is updated according to the First Law of Thermodynamics with entropy generation terms 
due to the presence of viscosity (see eq. EH1) : 



_n+l _ n 
'i+l/2 i + l/2 

St 



1 

2n 



(PI 



+1 

-1/2 



PI+1/2) 



V i+l/2 V i+l/2 

Si 



7) 



yn+1/2 .n+l/2, , 
V 1 + 1/2 J i+l 2 \ U >J> ' h 



(CIO) 
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where 



,>n-t-l/2 ( ■ , _ - 
fi+1/2 \V,J,r) = - 



Av. 



1+1/2 

f+1 



+ 1/2 
-1/2 



Ar™ 



' A(r? +1 )2 



+ 1/2 
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(C+l) 2 + (C) 2 



n / .n+1/2 

i+l/2_^ [ J i+ 1 



A 



■ n+l/2~ 
V i+1 



Ar n 
lAI i+1 



(C+l) 2 



(Cll) 



Note that the terms due to viscous dissipation are not time-centered. Equation (|C10|I may be solved for the updated 
value of the specific internal energy, s^f/g' ^ v nrs t substituting for the (as yet unknown) updated pressure, P"_^j 2 j 
according to the T-law EOS: 



P. 



1+1/2 



e 



i 

i+l/2/ v i+l/2- 



(C12) 



Finally, the pre ssure itself is updated using equation I|C12J| . In order to evolve the system from one timestep to the 
next, equations 1|C4|) - (|C12|) are used in the order as given above to updated each quantity. These equations only apply 
in the interior regions of the grid and some must be modified to accommodate boundary conditions on the cylinder 
axis (r = 0) and at the surface. These modifications result in only first order spatial accuracy near these boundaries. 

All of the calculations described in this paper were performed on single processors and none lasted more than a few 
hours. As an example, consider the rapidly rotating n — 3 cylinder with 77 = 0. This run required 2400 Lagrangian 
mass shells and each timestep required 7.295 ms of CPU time on a Pentium IV, 2.0 GHz processor. The total length 
of this run was 4.4 x 10 6 steps. Our runs required spatial resolution in the range of 400-2400 Lagrangian mass shells. 
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